Liana Merk
from bioscrape.types import Model
from bioscrape.simulator import py_simulate_model
from biocrnpyler import *
import numpy as np
import pylab as plt
import tqdm
import pandas as pd
import matplotlib as mpl
import holoviews as hv
import bokeh.plotting
import bokeh.io
from bokeh.models import Range1d
from bokeh.layouts import gridplot
bokeh.io.output_notebook()
hv.extension('bokeh')
from bokeh.models import PrintfTickFormatter
import seaborn as sns
sns.set()
%matplotlib inline
#Parameters": List of tuples [("paramter name" (string)": value (float))]
parameters = {"k1b": 100, "k1u": 1,
"k2b": 0.3, "k2u": 10,
"k3b": 0.3, "k3u": 10,
"k4b": 1.6, "k4u": 0.016,
"k5b": 0.0001, "k5u":1,
"k6b": 0.0001, "k6u": 0.1,
"k7b": 1, "k7u": 1,
"k8b": 0.0083, "k8u": 5,
"k9b": 0.3, "k9u": 3,
"k10b": 0.1, "k10u":1,
"k11b": 100, "k11u": 1,
"ktx": 0.09, "ktl": 0.05,
"k_dephos": 50, "delta": .9}
# Ribosome
R = Species("R")
# RNA Polymerase
P = Species("P")
# Consitutive Promoter P1d
P1d = Species("P1d")
# Polymerase bound to P1d
P_P1d = Species("P:P1d")
# P + P1d <-> P1d:P
R1 = Reaction([P, P1d], [P_P1d], k = parameters['k1b'], k_rev = parameters['k1u'])
# ttrS Transcript
ttrS_T = Species("ttrS_T")
# ttrR Transcript
ttrR_T = Species("ttrR_T")
# P1d:P -> ttrS_T + ttrR_T
Rtx1 = Reaction([P_P1d], [ttrS_T, ttrR_T, P, P1d], k = parameters['ktx'])
# ttrS trasncript bound to Ribosome
ttrS_T_R = Species("ttrS_T:R")
# ttrS
ttrS = Species("ttrS")
# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = parameters['k2u'])
# ttrS Translation
Rtl1 = Reaction([ttrS_T_R], [ttrS_T, ttrS, R], k = parameters['ktl'])
# ttrR trasncript bound to Ribosome
ttrR_T_R = Species("ttrR_T:R")
# ttrR
ttrR = Species("ttrR")
# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = parameters['k3u'])
# ttrR Translation
Rtl2 = Reaction([ttrR_T_R], [ttrR_T, ttrR, R], k = parameters['ktl'])
# Tetrathionate
tt = Species("tt")
# Phosphorylated ttrS
ttrS_P = Species("ttrS_P")
# ttrS + tt <-> ttrS_P + tt
R4 = Reaction([ttrS, tt], [ttrS_P, tt], k = parameters['k4b'], k_rev = parameters['k4u'])
# ttrR bound to ttrS
ttrR_ttrS = Species("ttrR_ttrS")
# ttrR + ttrS <-> ttrR:ttrS
R5 = Reaction([ttrS, ttrR], [ttrR_ttrS], k = parameters['k5b'], k_rev = parameters['k5u'])
# ttrR bound to phosphorylated ttrS
ttrR_ttrS_P = Species("ttrR_ttrS_P")
# ttrR + ttrS <-> ttrR:ttrS:P
R6 = Reaction([ttrS_P, ttrR], [ttrR_ttrS_P], k = parameters['k6b'], k_rev = parameters['k6u'])
# Phosphorylated ttrR
ttrR_P = Species("ttrR_P")
# ttrR:ttrS:P <-> ttrR:P + ttrS
R7 = Reaction([ttrR_ttrS_P], [ttrR_P, ttrS], k = parameters['k7b'], k_rev = parameters['k7u'])
# ttrR:P -> ttrR
rxn_dephos = Reaction([ttrR_P], [ttrR], k = parameters['k_dephos'])
# ttrR:P_dimer
ttrR_P_dimer = Species("ttrR:P_dimer")
# 2ttrR:P <-> ttrR:P_dimer
R8 = Reaction([ttrR_P, ttrR_P], [ttrR_P_dimer], k = parameters['k8b'], k_rev = parameters['k8u'])
# Inactive pttr
pttr0 = Species("pttr0")
# Inactive pttr
pttr1 = Species("pttr1")
# ttrR:P_dimer + pttr0 <-> pttr1
R9 = Reaction([ttrR_P_dimer, pttr0], [pttr1], k = parameters['k9b'], k_rev = parameters['k9u'])
# Active pttr to polymerase
pttr1_P = Species("pttr1_P")
# P + pttr1 <-> pttr1:P
Rtx2 = Reaction([P, pttr1], [pttr1_P], k = parameters['k10b'], k_rev = parameters['k10u'])
# Active pttr to polymerase
GFP_T = Species("GFP_T")
# P + pttr1 <-> pttr1:P
R10 = Reaction([pttr1_P], [GFP_T, P, pttr1], k = parameters['ktx'])
# Active pttr to polymerase
GFP_T_R = Species("GFP_T_R")
# GFP
GFP = Species("GFP")
# P + pttr1 <-> pttr1:P
R11 = Reaction([GFP_T, R], [GFP_T_R], k = parameters['k11b'], k_rev = parameters['k11u'])
# ttrS Translation
Rtl3 = Reaction([GFP_T_R], [GFP, GFP_T, R], k = parameters['ktl'])
# RNA deg
rxnd_ttrS_T = Reaction([ttrS_T], [], k = parameters['delta'])
rxnd_ttrR_T = Reaction([ttrR_T], [], k = parameters['delta'])
rxnd_GFP_T = Reaction([GFP_T], [], k = parameters['delta'])
# Protein Deg
rxnd_ttrS = Reaction([ttrS], [], k = parameters['delta'])
rxnd_ttrR = Reaction([ttrR], [], k = parameters['delta'])
rxnd_GFP = Reaction([GFP], [], k = parameters['delta'])
reactions = [R1, R2, R3, R4, R5, R6, R7, R8, R9, R10, R11,
Rtx1, Rtx2,
Rtl1, Rtl2, Rtl3,
rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
species = [R, P, P1d, P_P1d, ttrS_T, ttrR_T, ttrR_T_R, ttrS_T_R,
tt, ttrS_P, ttrR_ttrS, ttrR_ttrS_P, ttrR_P, ttrR_P_dimer,
pttr0, pttr1, pttr1_P, GFP_T, GFP_T_R, GFP]
# Create CRN
CRN_tt = ChemicalReactionNetwork(species = species, reactions = reactions)
print(CRN_tt.pretty_print(show_rates = True))
# An initial condition for each species (uninitialized species default to 0)
x0 = {
"tt":200,
"P1d":5,
"pttr0":5,
"P":25,
"R":100,
}
timepoints = np.linspace(0, 600, 100)
Results = CRN_tt.simulate_with_bioscrape(timepoints,x0)
Results.head()
def plot_holoviews_ttr_bioscrape(results_melted):
plot_group = []
for row in results_melted.iterrows():
sp = row[1]['species']
if sp in ["P1d", "P:P1d", "pttr0", "pttr1", "pttr1_P"]:
plot_group.append('promoters')
elif sp in ["ttrS_T", "ttrR_T", "ttrR_T:R", "ttrS_T:R", "GFP_T", "GFP_T_R"]:
plot_group.append('mRNA')
elif sp in ["ttrS_P", "ttrR_ttrS", "ttrR_ttrS_P", "ttrR_P", "ttrR:P_dimer"]:
plot_group.append('complexes')
elif sp in ["ttrR", "ttrS", "GFP"]:
plot_group.append('proteins')
elif sp in ["P", "R", "tt"]:
plot_group.append('machinery and signals')
else:
print(sp)
results_melted['plot_group'] = plot_group
plot_list = []
for group in results_melted['plot_group'].unique():
plot_list.append(hv.Curve(
data=results_melted.loc[results_melted['plot_group'] == group],
kdims=['time', 'value'],
vdims = ['species']
).groupby(
[ 'species']
).opts(
tools=['hover'],
# legend_position = 'right',
show_grid=True,
title = f'{group}',
width = 500,
#shared_axes=False
).overlay('species').opts(legend_position = 'right',))
return plot_list
results_melted = pd.melt(Results, var_name='species', id_vars=['time'])
results_melted.head()
plot_list = plot_holoviews_ttr_bioscrape(results_melted)
hv.Layout(plot_list).cols(2).opts(shared_axes=False)
def create_reactions(parameters):
# P + P1d <-> P1d:P
R1 = Reaction([P, P1d], [P_P1d], k = parameters['k1b'], k_rev = parameters['k1u'])
# P1d:P -> ttrS_T + ttrR_T
Rtx1 = Reaction([P_P1d], [ttrS_T, ttrR_T, P, P1d], k = parameters['ktx'])
# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = parameters['k2u'])
# ttrS Translation
Rtl1 = Reaction([ttrS_T_R], [ttrS_T, ttrS, R], k = parameters['ktl'])
# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = parameters['k3u'])
# ttrR Translation
Rtl2 = Reaction([ttrR_T_R], [ttrR_T, ttrR, R], k = parameters['ktl'])
# ttrS + tt <-> ttrS_P + tt
R4 = Reaction([ttrS, tt], [ttrS_P, tt], k = parameters['k4b'], k_rev = parameters['k4u'])
# ttrR + ttrS <-> ttrR:ttrS
R5 = Reaction([ttrS, ttrR], [ttrR_ttrS], k = parameters['k5b'], k_rev = parameters['k5u'])
# ttrR + ttrS <-> ttrR:ttrS:P
R6 = Reaction([ttrS_P, ttrR], [ttrR_ttrS_P], k = parameters['k6b'], k_rev = parameters['k6u'])
# ttrR:ttrS:P <-> ttrR:P + ttrS
R7 = Reaction([ttrR_ttrS_P], [ttrR_P, ttrS], k = parameters['k7b'], k_rev = parameters['k7u'])
# ttrR:P -> ttrR
rxn_dephos = Reaction([ttrR_P], [ttrR], k = parameters['k_dephos'])
# 2ttrR:P <-> ttrR:P_dimer
R8 = Reaction([ttrR_P, ttrR_P], [ttrR_P_dimer], k = parameters['k8b'], k_rev = parameters['k8u'])
# ttrR:P_dimer + pttr0 <-> pttr1
R9 = Reaction([ttrR_P_dimer, pttr0], [pttr1], k = parameters['k9b'], k_rev = parameters['k9u'])
# P + pttr1 <-> pttr1:P
Rtx2 = Reaction([P, pttr1], [pttr1_P], k = parameters['k10b'], k_rev = parameters['k10u'])
# P + pttr1 <-> pttr1:P
R10 = Reaction([pttr1_P], [GFP_T, P, pttr1], k = parameters['ktx'])
# P + pttr1 <-> pttr1:P
R11 = Reaction([GFP_T, R], [GFP_T_R], k = parameters['k11b'], k_rev = parameters['k11u'])
# ttrS Translation
Rtl3 = Reaction([GFP_T_R], [GFP, GFP_T, R], k = parameters['ktl'])
# RNA deg
rxnd_ttrS_T = Reaction([ttrS_T], [], k = parameters['delta'])
rxnd_ttrR_T = Reaction([ttrR_T], [], k = parameters['delta'])
rxnd_GFP_T = Reaction([GFP_T], [], k = parameters['delta'])
# Protein Deg
rxnd_ttrS = Reaction([ttrS], [], k = parameters['delta'])
rxnd_ttrR = Reaction([ttrR], [], k = parameters['delta'])
rxnd_GFP = Reaction([GFP], [], k = parameters['delta'])
reactions = [R1, R4, R5, R6, R7, R8, R9, R10, R11,
Rtx1, Rtx2,
Rtl1, Rtl2, Rtl3,
rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
return reactions
We will change k2u, and k3u, the reaction rates that define ribosome unbinding from ttrS and ttrR transcripts, respectively. We will save the simulations as csv's to plot in another notebook such that the styles match well.
all_sims = pd.DataFrame()
#Parameters": List of tuples [("paramter name" (string)": value (float))]
parameters = {"k1b": 100, "k1u": 1,
"k2b": 3, "k2u": 10,
"k3b": 0.00003, "k3u": 10,
"k4b": 5, "k4u": 0.5,
"k5b": 5, "k5u":6,
"k6b": 5, "k6u": 1,
"k7b": 0.01, "k7u": 1,
"k8b": 80, "k8u": 5,
"k9b": 0.004, "k9u": 3,
"k10b": 50, "k10u":1,
"k11b": 10, "k11u": 1,
"ktx": 0.1, "ktl": 1,
"k_dephos": 7, "delta": .9}
def plot_ttr_RBS_tuning(parameters, x0, k2u, k3u, title, timepoints):
experimental_ttr_concs = np.array([ 0, 0.1, 0.25, 0.5, 0.75, 1, 5, 10])
# Get color palette in order
r = list(bokeh.palettes.viridis(10))
r.reverse()
r = r[2:]
# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = k2u)
# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = k3u)
reactions = create_reactions(parameters)
reactions.append(R2)
reactions.append(R3)
species = [R, P, P1d, P_P1d, ttrS_T, ttrR_T, ttrR_T_R, ttrS_T_R,
tt, ttrS_P, ttrR_ttrS, ttrR_ttrS_P, ttrR_P, ttrR_P_dimer,
pttr0, pttr1, pttr1_P, GFP_T, GFP_T_R, GFP]
# Recompile the CRN with the right parameters
CRN_curr = ChemicalReactionNetwork(species = species, reactions = reactions)
p = bokeh.plotting.figure(width = 300, height= 250, x_axis_label='Time (hr)', y_axis_label = '[GFP]', title = title)
i = 0
for tt_conc in experimental_ttr_concs:
x0["tt"] = tt_conc
Results = CRN_curr.simulate_with_bioscrape(timepoints,x0)
all_sims[f'{title}_{tt_conc}'] = Results['GFP']
p.line(timepoints/60, Results['GFP'],color = r[i], legend_label = f'{tt_conc}', line_width = 2)
i+=1
p.legend.location = 'top_left'
if k2u != 20:
p.legend.visible = False
p.y_range = Range1d(-0.1,2.5)
return p
timepoints = np.linspace(0, 1400, 1400)
title_list = ['Very Low (LM15) Simulation', 'Low (LM14) Simulation', 'Medium (LM18) Simulation', 'High (LM19) Simulation']
x0 = {
"tt":200,
# PSC101 has copy number 5
"P1d":5,
"pttr0":5,
# Polymerase: 15 uM assuming ~10000 Polymerase molecules / E. coli with a volume 1 um^3
"P":15,
# Ribosome: 120 uM assuming ~72000 Ribosomes / E. coli with a volume 1 um^3
"R":120,
}
parameters = {"k1b": 10, "k1u": 0.0001,
"k2b": 0.3, "k2u": 10,
"k3b": 0.3, "k3u": 10,
"k4b": 1.6, "k4u": 0.016,
"k5b": 0.0001, "k5u":1,
"k6b": 0.0001, "k6u": 0.1,
"k7b": 1, "k7u": 1,
"k8b": 0.0083, "k8u": 0.5,
"k9b": 0.3, "k9u": 3,
"k10b": 10, "k10u":0.0001,
"k11b": 100, "k11u": 1,
"ktx": 0.1, "ktl":0.33,
"k_dephos": 50, "delta": .9}
# Define higher off rates
k2u = [20, 12, 11, 10]
k3u = [15, 8, 6, 5]
plot_list = []
for i in range(len(k3u)):
plot_list.append(plot_ttr_RBS_tuning(parameters, x0, k2u[i], k3u[i], title_list[i], timepoints))
bokeh.io.show(gridplot([plot_list]))
#Parameters": List of tuples [("paramter name" (string)": value (float))]
parameters = {"k1b": 100, "k1u": 1,
"k2b": 0.3, "k2u": 10,
"k3b": 0.3, "k3u": 10,
"k4b": 1.6, "k4u": 0.016,
"k5b": 0.0001, "k5u":1,
"k6b": 0.0001, "k6u": 0.1,
"k7b": 1, "k7u": 1,
"k8b": 0.0083, "k8u": 5,
"k9b": 0.3, "k9u": 3,
"k10b": 0.1, "k10u":1,
"k11b": 100, "k11u": 1,
"ktx": 0.09, "ktl": 0.05,
"k_dephos": 50, "delta": .9}
# Ribosome
R = Species("R")
# RNA Polymerase
P = Species("P")
# Consitutive Promoter P1d
P1d = Species("P1d")
# Polymerase bound to P1d
P_P1d = Species("P:P1d")
# P + P1d <-> P1d:P
R1 = Reaction([P, P1d], [P_P1d], k = parameters['k1b'], k_rev = parameters['k1u'])
# ttrS Transcript
ttrS_T = Species("ttrS_T")
# ttrR Transcript
ttrR_T = Species("ttrR_T")
# P1d:P -> ttrS_T + ttrR_T
Rtx1 = Reaction([P_P1d], [ttrS_T, ttrR_T, P, P1d], k = parameters['ktx'])
# ttrS trasncript bound to Ribosome
ttrS_T_R = Species("ttrS_T:R")
# ttrS
ttrS = Species("ttrS")
# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = parameters['k2u'])
# ttrS Translation
Rtl1 = Reaction([ttrS_T_R], [ttrS_T, ttrS, R], k = parameters['ktl'])
# ttrR trasncript bound to Ribosome
ttrR_T_R = Species("ttrR_T:R")
# ttrR
ttrR = Species("ttrR")
# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = parameters['k3u'])
# ttrR Translation
Rtl2 = Reaction([ttrR_T_R], [ttrR_T, ttrR, R], k = parameters['ktl'])
# Tetrathionate
tt = Species("tt")
# Phosphorylated ttrS
ttrS_P = Species("ttrS_P")
# ttrS + tt <-> ttrS_P + tt
R4 = Reaction([ttrS, tt], [ttrS_P, tt], k = parameters['k4b'], k_rev = parameters['k4u'])
# ttrR bound to ttrS
ttrR_ttrS = Species("ttrR_ttrS")
# ttrR + ttrS <-> ttrR:ttrS
R5 = Reaction([ttrS, ttrR], [ttrR_ttrS], k = parameters['k5b'], k_rev = parameters['k5u'])
# ttrR bound to phosphorylated ttrS
ttrR_ttrS_P = Species("ttrR_ttrS_P")
# ttrR + ttrS <-> ttrR:ttrS:P
R6 = Reaction([ttrS_P, ttrR], [ttrR_ttrS_P], k = parameters['k6b'], k_rev = parameters['k6u'])
# Phosphorylated ttrR
ttrR_P = Species("ttrR_P")
# ttrR:ttrS:P <-> ttrR:P + ttrS
R7 = Reaction([ttrR_ttrS_P], [ttrR_P, ttrS], k = parameters['k7b'], k_rev = parameters['k7u'])
# ttrR:P -> ttrR
rxn_dephos = Reaction([ttrR_P], [ttrR], k = parameters['k_dephos'])
# ttrR:P_dimer
ttrR_P_dimer = Species("ttrR:P_dimer")
# 2ttrR:P <-> ttrR:P_dimer
R8 = Reaction([ttrR_P, ttrR_P], [ttrR_P_dimer], k = parameters['k8b'], k_rev = parameters['k8u'])
# Inactive pttr
pttr0 = Species("pttr0")
# Inactive pttr
pttr1 = Species("pttr1")
# ttrR:P_dimer + pttr0 <-> pttr1
R9 = Reaction([ttrR_P_dimer, pttr0], [pttr1], k = parameters['k9b'], k_rev = parameters['k9u'])
# Active pttr to polymerase
pttr1_P = Species("pttr1_P")
# P + pttr1 <-> pttr1:P
Rtx2 = Reaction([P, pttr1], [pttr1_P], k = parameters['k10b'], k_rev = parameters['k10u'])
# Active pttr to polymerase
GFP_T = Species("GFP_T")
# P + pttr1 <-> pttr1:P
R10 = Reaction([pttr1_P], [GFP_T, P, pttr1], k = parameters['ktx'])
# Active pttr to polymerase
GFP_T_R = Species("GFP_T_R")
# GFP
GFP = Species("GFP")
# P + pttr1 <-> pttr1:P
R11 = Reaction([GFP_T, R], [GFP_T_R], k = parameters['k11b'], k_rev = parameters['k11u'])
# ttrS Translation
Rtl3 = Reaction([GFP_T_R], [GFP, GFP_T, R], k = parameters['ktl'])
# RNA deg
rxnd_ttrS_T = Reaction([ttrS_T], [], k = parameters['delta'])
rxnd_ttrR_T = Reaction([ttrR_T], [], k = parameters['delta'])
rxnd_GFP_T = Reaction([GFP_T], [], k = parameters['delta'])
# Protein Deg
rxnd_ttrS = Reaction([ttrS], [], k = parameters['delta'])
rxnd_ttrR = Reaction([ttrR], [], k = parameters['delta'])
rxnd_GFP = Reaction([GFP], [], k = parameters['delta'])
reactions = [R1, R2, R3, R4, R5, R6, R7, R8, R9, R10, R11,
Rtx1, Rtx2,
Rtl1, Rtl2, Rtl3,
rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
species = [R, P, P1d, P_P1d, ttrS_T, ttrR_T, ttrR_T_R, ttrS_T_R,
tt, ttrS_P, ttrR_ttrS, ttrR_ttrS_P, ttrR_P, ttrR_P_dimer,
pttr0, pttr1, pttr1_P, GFP_T, GFP_T_R, GFP]
# Create CRN
CRN_tt = ChemicalReactionNetwork(species = species, reactions = reactions)
print(CRN_tt.pretty_print(show_rates = True))
# An initial condition for each species (uninitialized species default to 0)
x0 = {
"tt":200,
"P1d":5,
"pttr0":5,
"P":25,
"R":100,
}
timepoints = np.linspace(0, 600, 100)
Results = CRN_tt.simulate_with_bioscrape(timepoints,x0)
Results.head()
def plot_holoviews_ttr_bioscrape(results_melted):
plot_group = []
for row in results_melted.iterrows():
sp = row[1]['species']
if sp in ["P1d", "P:P1d", "pttr0", "pttr1", "pttr1_P"]:
plot_group.append('promoters')
elif sp in ["ttrS_T", "ttrR_T", "ttrR_T:R", "ttrS_T:R", "GFP_T", "GFP_T_R"]:
plot_group.append('mRNA')
elif sp in ["ttrS_P", "ttrR_ttrS", "ttrR_ttrS_P", "ttrR_P", "ttrR:P_dimer"]:
plot_group.append('complexes')
elif sp in ["ttrR", "ttrS", "GFP"]:
plot_group.append('proteins')
elif sp in ["P", "R", "tt"]:
plot_group.append('machinery and signals')
else:
print(sp)
results_melted['plot_group'] = plot_group
plot_list = []
for group in results_melted['plot_group'].unique():
plot_list.append(hv.Curve(
data=results_melted.loc[results_melted['plot_group'] == group],
kdims=['time', 'value'],
vdims = ['species']
).groupby(
[ 'species']
).opts(
tools=['hover'],
# legend_position = 'right',
show_grid=True,
title = f'{group}',
width = 500,
#shared_axes=False
).overlay('species').opts(legend_position = 'right',))
return plot_list
results_melted = pd.melt(Results, var_name='species', id_vars=['time'])
results_melted.head()
plot_list = plot_holoviews_ttr_bioscrape(results_melted)
hv.Layout(plot_list).cols(2).opts(shared_axes=False)
def create_reactions(parameters):
# P + P1d <-> P1d:P
R1 = Reaction([P, P1d], [P_P1d], k = parameters['k1b'], k_rev = parameters['k1u'])
# P1d:P -> ttrS_T + ttrR_T
Rtx1 = Reaction([P_P1d], [ttrS_T, ttrR_T, P, P1d], k = parameters['ktx'])
# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = parameters['k2u'])
# ttrS Translation
Rtl1 = Reaction([ttrS_T_R], [ttrS_T, ttrS, R], k = parameters['ktl'])
# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = parameters['k3u'])
# ttrR Translation
Rtl2 = Reaction([ttrR_T_R], [ttrR_T, ttrR, R], k = parameters['ktl'])
# ttrS + tt <-> ttrS_P + tt
R4 = Reaction([ttrS, tt], [ttrS_P, tt], k = parameters['k4b'], k_rev = parameters['k4u'])
# ttrR + ttrS <-> ttrR:ttrS
R5 = Reaction([ttrS, ttrR], [ttrR_ttrS], k = parameters['k5b'], k_rev = parameters['k5u'])
# ttrR + ttrS <-> ttrR:ttrS:P
R6 = Reaction([ttrS_P, ttrR], [ttrR_ttrS_P], k = parameters['k6b'], k_rev = parameters['k6u'])
# ttrR:ttrS:P <-> ttrR:P + ttrS
R7 = Reaction([ttrR_ttrS_P], [ttrR_P, ttrS], k = parameters['k7b'], k_rev = parameters['k7u'])
# ttrR:P -> ttrR
rxn_dephos = Reaction([ttrR_P], [ttrR], k = parameters['k_dephos'])
# 2ttrR:P <-> ttrR:P_dimer
R8 = Reaction([ttrR_P, ttrR_P], [ttrR_P_dimer], k = parameters['k8b'], k_rev = parameters['k8u'])
# ttrR:P_dimer + pttr0 <-> pttr1
R9 = Reaction([ttrR_P_dimer, pttr0], [pttr1], k = parameters['k9b'], k_rev = parameters['k9u'])
# P + pttr1 <-> pttr1:P
Rtx2 = Reaction([P, pttr1], [pttr1_P], k = parameters['k10b'], k_rev = parameters['k10u'])
# P + pttr1 <-> pttr1:P
R10 = Reaction([pttr1_P], [GFP_T, P, pttr1], k = parameters['ktx'])
# P + pttr1 <-> pttr1:P
R11 = Reaction([GFP_T, R], [GFP_T_R], k = parameters['k11b'], k_rev = parameters['k11u'])
# ttrS Translation
Rtl3 = Reaction([GFP_T_R], [GFP, GFP_T, R], k = parameters['ktl'])
# RNA deg
rxnd_ttrS_T = Reaction([ttrS_T], [], k = parameters['delta'])
rxnd_ttrR_T = Reaction([ttrR_T], [], k = parameters['delta'])
rxnd_GFP_T = Reaction([GFP_T], [], k = parameters['delta'])
# Protein Deg
rxnd_ttrS = Reaction([ttrS], [], k = parameters['delta'])
rxnd_ttrR = Reaction([ttrR], [], k = parameters['delta'])
rxnd_GFP = Reaction([GFP], [], k = parameters['delta'])
reactions = [R1, R4, R5, R6, R7, R8, R9, R10, R11,
Rtx1, Rtx2,
Rtl1, Rtl2, Rtl3,
rxnd_ttrS_T, rxnd_ttrR_T, rxnd_GFP_T,
rxnd_ttrS, rxnd_ttrR, rxnd_GFP]
return reactions
We will change k2u, and k3u, the reaction rates that define ribosome unbinding from ttrS and ttrR transcripts, respectively. We will save the simulations as csv's to plot in another notebook such that the styles match well.
all_sims = pd.DataFrame()
#Parameters": List of tuples [("paramter name" (string)": value (float))]
parameters = {"k1b": 100, "k1u": 1,
"k2b": 3, "k2u": 10,
"k3b": 0.00003, "k3u": 10,
"k4b": 5, "k4u": 0.5,
"k5b": 5, "k5u":6,
"k6b": 5, "k6u": 1,
"k7b": 0.01, "k7u": 1,
"k8b": 80, "k8u": 5,
"k9b": 0.004, "k9u": 3,
"k10b": 50, "k10u":1,
"k11b": 10, "k11u": 1,
"ktx": 0.1, "ktl": 1,
"k_dephos": 7, "delta": .9}
def plot_ttr_RBS_tuning(parameters, x0, k2u, k3u, title, timepoints):
experimental_ttr_concs = np.array([ 0, 0.1, 0.25, 0.5, 0.75, 1, 5, 10])
# Get color palette in order
r = list(bokeh.palettes.viridis(10))
r.reverse()
r = r[2:]
# ttrS_T + R <-> ttrS_T:R
R2 = Reaction([ttrS_T, R], [ttrS_T_R], k = parameters['k2b'], k_rev = k2u)
# ttrR_T + R <-> ttrR_T:R
R3 = Reaction([ttrR_T, R], [ttrR_T_R], k = parameters['k3b'], k_rev = k3u)
reactions = create_reactions(parameters)
reactions.append(R2)
reactions.append(R3)
species = [R, P, P1d, P_P1d, ttrS_T, ttrR_T, ttrR_T_R, ttrS_T_R,
tt, ttrS_P, ttrR_ttrS, ttrR_ttrS_P, ttrR_P, ttrR_P_dimer,
pttr0, pttr1, pttr1_P, GFP_T, GFP_T_R, GFP]
# Recompile the CRN with the right parameters
CRN_curr = ChemicalReactionNetwork(species = species, reactions = reactions)
p = bokeh.plotting.figure(width = 300, height= 250, x_axis_label='Time (hr)', y_axis_label = '[GFP]', title = title)
i = 0
for tt_conc in experimental_ttr_concs:
x0["tt"] = tt_conc
Results = CRN_curr.simulate_with_bioscrape(timepoints,x0)
all_sims[f'{title}_{tt_conc}'] = Results['GFP']
p.line(timepoints/60, Results['GFP'],color = r[i], legend_label = f'{tt_conc}', line_width = 2)
i+=1
p.legend.location = 'top_left'
if k2u != 20:
p.legend.visible = False
p.y_range = Range1d(-0.1,2.5)
return p
timepoints = np.linspace(0, 1400, 1400)
title_list = ['Very Low (LM15) Simulation', 'Low (LM14) Simulation', 'Medium (LM18) Simulation', 'High (LM19) Simulation']
x0 = {
"tt":200,
# PSC101 has copy number 5
"P1d":5,
"pttr0":5,
# Polymerase: 15 uM assuming ~10000 Polymerase molecules / E. coli with a volume 1 um^3
"P":15,
# Ribosome: 120 uM assuming ~72000 Ribosomes / E. coli with a volume 1 um^3
"R":120,
}
parameters = {"k1b": 10, "k1u": 0.0001,
"k2b": 0.3, "k2u": 10,
"k3b": 0.3, "k3u": 10,
"k4b": 1.6, "k4u": 0.016,
"k5b": 0.0001, "k5u":1,
"k6b": 0.0001, "k6u": 0.1,
"k7b": 1, "k7u": 1,
"k8b": 0.0083, "k8u": 0.5,
"k9b": 0.3, "k9u": 3,
"k10b": 10, "k10u":0.0001,
"k11b": 100, "k11u": 1,
"ktx": 0.1, "ktl":0.33,
"k_dephos": 50, "delta": .9}
# Define higher off rates
k2u = [20, 12, 11, 10]
k3u = [15, 8, 6, 5]
plot_list = []
for i in range(len(k3u)):
plot_list.append(plot_ttr_RBS_tuning(parameters, x0, k2u[i], k3u[i], title_list[i], timepoints))
bokeh.io.show(gridplot([plot_list]))
all_sims['Time (hr)'] = timepoints
all_sims.to_csv('./data/all_sims.csv', index = False)
This $\sigma 54$-dependent split activator requires two proteins, hrpR and hrpS, to function. For now, let's just spike in the two proteins and see how our and gate responds.
txtl = CRNLab("GFP")
txtl.mixture("mixture1", extract = "TxTlExtract", mixture_volume = 1e-6, mixture_parameters = './parameters/BasicExtract.tsv')
#Define the activator and repressor species as proteins
prot1 = Protein("hrpR")
prot2 = Protein("hrpS")
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
#TxTl Extract is a Mixture with more complex internal models
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
Re1.head()
Now, we will make sure the AND gate (combinatorial promoter) is working by varying protein concentrations.
Let's create a plotting function. We will plot a heatmap of protein GFP at the last slice of the simulation. Because we have no degradation of protein, the last timepoint will be the max. We will be able to pass in normalize, which divides all values by the maximum of all simulated values, or max_to_1, which will set the scale bar to be 0 to 1. This basically causes plots that might've been
def simulate_heatmap(inducer_1_name,
inducer_1_values,
inducer_2_name,
inducer_2_values,
title,
x0=x0,
normalize=True,
max_to_1 = False,
CRN=CRN_extract_1,
timepoints=timepoints):
GFP_max = pd.DataFrame()
for conc_1 in inducer_1_values:
x0[inducer_1_name] = conc_1
for conc_2 in inducer_2_values:
x0[inducer_2_name] = conc_2 #Change my initial condition dictionary
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
GFP_max = GFP_max.append({inducer_1_name:conc_1,
inducer_2_name:conc_2,
'GFP_max': Re1["protein_GFP"].values[-1]}, ignore_index=True)
data = pd.pivot_table(data = GFP_max, index = inducer_1_name,
columns = inducer_2_name,
values = 'GFP_max')
fig, ax = plt.subplots(figsize=(6, 4.5))
if normalize:
data /= np.max(GFP_max['GFP_max'])
if max_to_1:
heat_map = sns.heatmap(data,
cmap = 'viridis',
linewidths=.1, vmin = 0, vmax = 1,
cbar_kws={'label': 'Normalized GFP'})
else:
heat_map = sns.heatmap(data,
cmap = 'viridis',
linewidths=.2,
cbar_kws={'label': 'Normalized GFP'})
sns.set(font_scale=1.4)
# plt.tick_params(axis='both', labelsize=20)
plt.xticks(fontsize=18, fontstyle='normal', rotation = 90)
plt.yticks(fontsize=18, fontstyle='normal')
plt.xlabel(inducer_1_name, fontsize=18, fontweight='bold')
plt.ylabel(inducer_2_name, fontsize=18, fontweight='bold')
heat_map.set_title(title)
fig.tight_layout()
return fig, heat_map
Let's say hrpR mutated and now it can recruit RNAP better than it did beter. We will change phrpL_hrpR_RNAP ku = 5 instead of 500.
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and_leak5_hrpR.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 1x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)
heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_1.png", bbox_inches = "tight")
plt.show()
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and_leak50_hrpR.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 10x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)
heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_2.png", bbox_inches = "tight")
plt.show()
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 100x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)
heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_3.png", bbox_inches = "tight")
plt.show()
class LacIPromoter2(Promoter):
def __init__(self, name = "pLac", IPTG = "IPTG",
lacI = "lacI", leak = True, assembly = None,
transcript = None, length = 0, mechanisms = {},
parameters = {}, **keywords):
Promoter.__init__(self, name=name, transcript=transcript, assembly=assembly, length=length, parameters=parameters, **keywords)
#RegulatedPromoter.__init__(self = self, name = name, regulators = regulators)
# Transcription is already included in promoter class
self.default_mechanisms = {"binding": One_Step_Cooperative_Binding()}
self.leak = leak
self.IPTG = self.set_species(IPTG, material_type = "inducer")
self.lacI = self.set_species(lacI, material_type = "protein")
self.plac_2_lacI = None
self.IPTG_2_lacI = None
def update_species(self):
mech_tx = self.mechanisms['transcription']
mech_b = self.mechanisms['binding']
species = []
self.complexes = []
"""A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
2lacI <--> 2xlacI
2xlacI + pLac <--> 2xlacI:pLac
"""
species_b = mech_b.update_species(self.lacI, self.assembly.dna, n_mer_species=self.lacI, component = self, part_id = self.assembly.dna.name+'_2x_'+self.lacI.name, cooperativity=2)
species += species_b
#print('plac stuff \n',species_b)
for s in species_b:
if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.assembly.dna in s.species:
self.plac_2_lacI = s
if self.leak is not False:
species += mech_tx.update_species(dna = self.plac_2_lacI, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
# IPTG can bind to lacI to sequester it
# give citation for 1 being enough
"""A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
2lacI <--> 2xlacI
2xlacI + IPTG <--> 2xlacI:IPTG
"""
species_b = mech_b.update_species(self.lacI, self.IPTG, component = self, part_id = self.IPTG.name+'_2x_'+self.lacI.name, cooperativity=2)
#print('iptg stuff \n', species_b)
species += species_b
for s in species_b:
if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.IPTG in s.species:
self.IPTG_2_lacI = s
# Unbound transcription
species += mech_tx.update_species(dna = self.assembly.dna, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
return species
def update_reactions(self):
print('r')
mech_tx = self.mechanisms['transcription']
mech_b = self.mechanisms['binding']
if self.IPTG_2_lacI == None or self.plac_2_lacI == None:
self.update_species()
if self.leak != False:
# Will need transcription, plac_2x_lacI, ktx, ku, kb
reactions += mech_tx.update_reactions(dna = self.plac_2_lacI, component = self, part_id = self.plac_2_lacI.name, \
transcript = self.transcript, protein = self.assembly.protein)
reactions = []
# Repression
# Will need two_step_cooperative_binding, plac_2x_lacI, kb1 ku1, kb2, ku2
reactions += mech_b.update_reactions(self.lacI, self.assembly.dna, component = self, \
part_id = self.plac_2_lacI.name, cooperativity = 2)
# LacI binding to IPTG
# Will need two_step_cooperative_binding, IPTG_2x_lacI, kb1 ku1, kb2, ku2
# .name gets rid of the complex
# .type will give you complex, species, etc
# .attributes gives you properties like phosphorylation state
reactions += mech_b.update_reactions(self.lacI, self.IPTG, component = self, part_id = self.IPTG_2_lacI.name, cooperativity=2)
'''
Any bit of IPTG destroys the bound structure
IPTG + plac_2xlacI --> IPTG:2xlacI + plac
'''
if isinstance(mech_b, Two_Step_Cooperative_Binding):
kb = self.get_parameter("kb2", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
else:
kb = self.get_parameter("kb", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
print(self.IPTG_2_lacI.name, kb)
reactions.append(Reaction([self.IPTG, self.plac_2_lacI], [self.IPTG_2_lacI, self.assembly.dna], k=kb))
# Unbound transcription plac, kb, ku, ktx
reactions += mech_tx.update_reactions(dna = self.assembly.dna, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
return reactions
def param_printer2(x0, parameters):
inits = pd.DataFrame(list(x0.items()),columns = ['param','value'])
df = pd.DataFrame(list(parameters.items()),columns = ['param','value'])
col_part_id = []
mech = []
for val in df['param']:
if isinstance(val, tuple):
col_part_id.append(val[0])
mech.append(val[1])
else:
col_part_id.append(' ')
mech.append(val)
df['param'] = mech
df['part'] = col_part_id
df.reindex(['part','param','value'],axis=1)
df[''] = ['|']*len(col_part_id)
return pd.concat([df, inits], axis=1).fillna("")
parameters = {
('pttr_phosphate_P_protein_ttrR', 'cooperativity'): 2.0,
('pttr_phosphate_P_protein_ttrR', 'ktx'): 0.17025, #this is the transcription rate
('pttr_phosphate_P_protein_ttrR', 'ku'): 11.01321586, #this is the polymerase unbinding rate
('transcription_mm', 'ktx'): .01, #These are the parameters for transcription leak
('transcription_mm','kb'):100, #default binding rate!
#('ttrS-P',"kdephos"):1000, #auto dephosphorylation
('ligand_tt_protein_ttrS-P','kdephos'):1000, #dephosphorylation rate of ttrS
#('ttrS-P',"kphos"):10, #ttrs phosphorylates itself automatically. This is a proxy for sensing
('ligand_tt_protein_ttrS-P','kphos'):10, #equivalent to kphos, from before
#('ttrS-ttrR','ktransfer'):5, #rate of transfer from ttrS to ttrR
('ligand_tt_protein_ttrS-ttrR', 'ktransfer'):5, #this is the new transfer rate
('ttrR-P',"kdephos"):1000, #auto dephosphorylation
"kb":100, "ku":10, "ktx":.05, "ktl":.2, "kdeg":3, 'cooperativity':2, #some default parameters
('phrpL_hrpR', 'cooperativity'): 1.0, #default is 2, but these proteins bind at cooperativity = 1
('phrpL_hrpS', 'cooperativity'): 1.0,
('phrpL_hrpR_hrpS_RNAP', 'ktx'): 0.17025, #this is the transcription rate
('phrpL_hrpR_hrpS_RNAP', 'ku'): 11.01321586, #this is the polymerase unbinding rate
('translation_mm', 'B0030', 'ku'): 888, #Unbinding
('translation_mm', 'B0030', 'ktl'): .15, #Translation Rate
('transcription_mm', 'ktx'): .01, #These are the parameters for transcription leak
('transcription_mm', 'ku'): 100,
('phrpL_hrpR','ku'):100,
('phrpL_hrpS','ku'):100,
('phrpL_hrpS_RNAP', 'ku'): 1000.,#we'll say that leak has 100x the unbinding rate
('phrpL_hrpR_RNAP', 'ku'): 1000,
('phrpL','ku'):100000,
('binding', 'pLac_lacI','ktx'):0.001,
('plac_lacI','kb'):10,
('plac_lacI','ku'):0.001,
('binding', 'plac_hrpS_2x_lacI','kb'):8888,
("simple_transcription", 'plac','ktx'):0.1,
#Make a lot of repressor
('translation_mm', 'rbsL', 'ku'): 0.0001, #Unbinding
('translation_mm', 'rbsL', 'ktl'): 2, #Translation Rate
('plac_hrpS','kb'):100,
('plac_hrpS','ku'):0.0001
}
#P_rbsS_ttrS
constit_ttrS = DNAassembly(name = "constit_ttrS", promoter = "P", rbs = "rbsS", protein = "ttrS")
#P_rbsR_ttrR
constit_ttrR = DNAassembly(name = "constit_ttrR", promoter = "P", rbs = "rbsR", protein = "ttrR")
#P_arl_lacI
constit_lacI = DNAassembly(name = "constit_lacI", promoter = "P", rbs = "rbsL", protein = "lacI")
#P_arl_lacI
constit_mScarlet = DNAassembly(name = "constit_mScarlet", promoter = "P", rbs = "B0034", protein = "mScarlet")
# plac_hrpS
IPTG = Species("IPTG", material_type='inducer')
#IPTG = Protein("IPTG")
lac_p = LacIPromoter2(name="pLac", IPTG = "IPTG", lacI = "lacI", leak=False)
plac_hrpS = DNAassembly(name = "plac_hrpS", promoter=lac_p, rbs="rbsH",protein="hrpS")
# pttr_hrpR
tt = Species("tt",material_type="ligand")
ttrSbound = ChemicalComplex([Protein("ttrS"),tt])
ttrS3 = PhosphoTransferase(ttrSbound.species,targets=["ttrR"])
ttrR = PhosphoProtein("ttrR")
# Need to change from UTR1
pttr_hrpR = DNAassembly(name = "pttr_hrpR",promoter=RegulatedPromoter("pttr",[ttrR.get_phosphorylated_species()],leak=False),rbs="UTR1",protein="hrpR")
hrpL_gfp = DNAassembly("hrpL_gfp",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=False),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli", components = [pttr_hrpR,ttrSbound,ttrS3,ttrR,constit_ttrS,constit_ttrR,
constit_lacI, plac_hrpS, hrpL_gfp, constit_mScarlet],parameters=parameters)
CRN_extract_1 = extract_1_TXTL.compile_crn()
param_printer2(x0, parameters)
x0 = {"dna_pttr_hrpR":1,
"dna_constit_ttrS":1,
"dna_constit_ttrR":1,
"dna_constit_lacI":1,
"dna_plac_hrpS":1,
"dna_hrpL_gfp":1,
"dna_constit_mScarlet":5,
"phosphate_P":100,
"protein_RNAP":15,
"protein_RNAase":30,
"protein_Ribo":120,
"protein_lacI": 5}
IPTG_um_list = np.logspace(-3, 2, 6).round(decimals=6)
tt_um_list = np.logspace(-3, 2, 6).round(decimals=6)
fig, heat_map = simulate_heatmap('ligand_tt',tt_um_list,'inducer_IPTG',IPTG_um_list,x0=x0,normalize=True,title='Compiled Full AND gate \n lacI Spike', CRN=CRN_extract_1)
heat_map.set(xlabel='IPTG (uM)', ylabel='Tetrathionate (uM)')
heat_map.invert_yaxis()
plt.yticks(fontsize=18, fontstyle='normal', rotation = 0)
fig.savefig("./figures/fig6/fig6d.png")
plt.show()
all_sims['Time (hr)'] = timepoints
all_sims.to_csv('./data/all_sims.csv', index = False)
This $\sigma 54$-dependent split activator requires two proteins, hrpR and hrpS, to function. For now, let's just spike in the two proteins and see how our and gate responds.
txtl = CRNLab("GFP")
txtl.mixture("mixture1", extract = "TxTlExtract", mixture_volume = 1e-6, mixture_parameters = './parameters/BasicExtract.tsv')
#Define the activator and repressor species as proteins
prot1 = Protein("hrpR")
prot2 = Protein("hrpS")
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
#TxTl Extract is a Mixture with more complex internal models
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
Re1.head()
Now, we will make sure the AND gate (combinatorial promoter) is working by varying protein concentrations.
Let's create a plotting function. We will plot a heatmap of protein GFP at the last slice of the simulation. Because we have no degradation of protein, the last timepoint will be the max. We will be able to pass in normalize, which divides all values by the maximum of all simulated values, or max_to_1, which will set the scale bar to be 0 to 1. This basically causes plots that might've been
def simulate_heatmap(inducer_1_name,
inducer_1_values,
inducer_2_name,
inducer_2_values,
title,
x0=x0,
normalize=True,
max_to_1 = False,
CRN=CRN_extract_1,
timepoints=timepoints):
GFP_max = pd.DataFrame()
for conc_1 in inducer_1_values:
x0[inducer_1_name] = conc_1
for conc_2 in inducer_2_values:
x0[inducer_2_name] = conc_2 #Change my initial condition dictionary
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
GFP_max = GFP_max.append({inducer_1_name:conc_1,
inducer_2_name:conc_2,
'GFP_max': Re1["protein_GFP"].values[-1]}, ignore_index=True)
data = pd.pivot_table(data = GFP_max, index = inducer_1_name,
columns = inducer_2_name,
values = 'GFP_max')
fig, ax = plt.subplots(figsize=(6, 4.5))
if normalize:
data /= np.max(GFP_max['GFP_max'])
if max_to_1:
heat_map = sns.heatmap(data,
cmap = 'viridis',
linewidths=.1, vmin = 0, vmax = 1,
cbar_kws={'label': 'Normalized GFP'})
else:
heat_map = sns.heatmap(data,
cmap = 'viridis',
linewidths=.2,
cbar_kws={'label': 'Normalized GFP'})
sns.set(font_scale=1.4)
# plt.tick_params(axis='both', labelsize=20)
plt.xticks(fontsize=18, fontstyle='normal', rotation = 90)
plt.yticks(fontsize=18, fontstyle='normal')
plt.xlabel(inducer_1_name, fontsize=18, fontweight='bold')
plt.ylabel(inducer_2_name, fontsize=18, fontweight='bold')
heat_map.set_title(title)
fig.tight_layout()
return fig, heat_map
Let's say hrpR mutated and now it can recruit RNAP better than it did beter. We will change phrpL_hrpR_RNAP ku = 5 instead of 500.
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and_leak5_hrpR.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 1x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)
heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_1.png", bbox_inches = "tight")
plt.show()
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and_leak50_hrpR.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 10x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)
heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_2.png", bbox_inches = "tight")
plt.show()
hrpS_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
hrpR_values = [0.0001, 0.001, 0.01, 0.1, 0.5, 1]
dna = DNAassembly("comb_promoter",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=True),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli extract 1", components = [dna,prot1,prot2], parameter_file = "./parameters/parameters_and.txt")
CRN_extract_1 = extract_1_TXTL.compile_crn()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_comb_promoter":5.0, "protein_hrpS":100, "protein_hrpR":100}
fig, heat_map = simulate_heatmap('protein_hrpS',hrpS_values,'protein_hrpR',hrpR_values,x0=x0,normalize=True,title='Simulated Protein-Based \n AND gate 100x $k_u^{hrpR}$', CRN=CRN_extract_1)
heat_map.figure.axes[-1].yaxis.label.set_size(20)
heat_map.set(xlabel='hrpS protein', ylabel='hrpR protein')
heat_map.invert_yaxis()
fig.savefig("./figures/fig5/fig5c_3.png", bbox_inches = "tight")
plt.show()
class LacIPromoter2(Promoter):
def __init__(self, name = "pLac", IPTG = "IPTG",
lacI = "lacI", leak = True, assembly = None,
transcript = None, length = 0, mechanisms = {},
parameters = {}, **keywords):
Promoter.__init__(self, name=name, transcript=transcript, assembly=assembly, length=length, parameters=parameters, **keywords)
#RegulatedPromoter.__init__(self = self, name = name, regulators = regulators)
# Transcription is already included in promoter class
self.default_mechanisms = {"binding": One_Step_Cooperative_Binding()}
self.leak = leak
self.IPTG = self.set_species(IPTG, material_type = "inducer")
self.lacI = self.set_species(lacI, material_type = "protein")
self.plac_2_lacI = None
self.IPTG_2_lacI = None
def update_species(self):
mech_tx = self.mechanisms['transcription']
mech_b = self.mechanisms['binding']
species = []
self.complexes = []
"""A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
2lacI <--> 2xlacI
2xlacI + pLac <--> 2xlacI:pLac
"""
species_b = mech_b.update_species(self.lacI, self.assembly.dna, n_mer_species=self.lacI, component = self, part_id = self.assembly.dna.name+'_2x_'+self.lacI.name, cooperativity=2)
species += species_b
#print('plac stuff \n',species_b)
for s in species_b:
if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.assembly.dna in s.species:
self.plac_2_lacI = s
if self.leak is not False:
species += mech_tx.update_species(dna = self.plac_2_lacI, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
# IPTG can bind to lacI to sequester it
# give citation for 1 being enough
"""A reaction where n binders (s1) bind to 1 bindee (s2) in two steps
2lacI <--> 2xlacI
2xlacI + IPTG <--> 2xlacI:IPTG
"""
species_b = mech_b.update_species(self.lacI, self.IPTG, component = self, part_id = self.IPTG.name+'_2x_'+self.lacI.name, cooperativity=2)
#print('iptg stuff \n', species_b)
species += species_b
for s in species_b:
if isinstance(s, ComplexSpecies) and s.species.count(self.lacI) == 2 and self.IPTG in s.species:
self.IPTG_2_lacI = s
# Unbound transcription
species += mech_tx.update_species(dna = self.assembly.dna, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
return species
def update_reactions(self):
print('r')
mech_tx = self.mechanisms['transcription']
mech_b = self.mechanisms['binding']
if self.IPTG_2_lacI == None or self.plac_2_lacI == None:
self.update_species()
if self.leak != False:
# Will need transcription, plac_2x_lacI, ktx, ku, kb
reactions += mech_tx.update_reactions(dna = self.plac_2_lacI, component = self, part_id = self.plac_2_lacI.name, \
transcript = self.transcript, protein = self.assembly.protein)
reactions = []
# Repression
# Will need two_step_cooperative_binding, plac_2x_lacI, kb1 ku1, kb2, ku2
reactions += mech_b.update_reactions(self.lacI, self.assembly.dna, component = self, \
part_id = self.plac_2_lacI.name, cooperativity = 2)
# LacI binding to IPTG
# Will need two_step_cooperative_binding, IPTG_2x_lacI, kb1 ku1, kb2, ku2
# .name gets rid of the complex
# .type will give you complex, species, etc
# .attributes gives you properties like phosphorylation state
reactions += mech_b.update_reactions(self.lacI, self.IPTG, component = self, part_id = self.IPTG_2_lacI.name, cooperativity=2)
'''
Any bit of IPTG destroys the bound structure
IPTG + plac_2xlacI --> IPTG:2xlacI + plac
'''
if isinstance(mech_b, Two_Step_Cooperative_Binding):
kb = self.get_parameter("kb2", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
else:
kb = self.get_parameter("kb", part_id = self.IPTG_2_lacI.name, mechanism = mech_b)
print(self.IPTG_2_lacI.name, kb)
reactions.append(Reaction([self.IPTG, self.plac_2_lacI], [self.IPTG_2_lacI, self.assembly.dna], k=kb))
# Unbound transcription plac, kb, ku, ktx
reactions += mech_tx.update_reactions(dna = self.assembly.dna, component = self, part_id = self.name, \
transcript = self.transcript, protein = self.assembly.protein)
return reactions
def param_printer2(x0, parameters):
inits = pd.DataFrame(list(x0.items()),columns = ['param','value'])
df = pd.DataFrame(list(parameters.items()),columns = ['param','value'])
col_part_id = []
mech = []
for val in df['param']:
if isinstance(val, tuple):
col_part_id.append(val[0])
mech.append(val[1])
else:
col_part_id.append(' ')
mech.append(val)
df['param'] = mech
df['part'] = col_part_id
df.reindex(['part','param','value'],axis=1)
df[''] = ['|']*len(col_part_id)
return pd.concat([df, inits], axis=1).fillna("")
parameters = {
('pttr_phosphate_P_protein_ttrR', 'cooperativity'): 2.0,
('pttr_phosphate_P_protein_ttrR', 'ktx'): 0.17025, #this is the transcription rate
('pttr_phosphate_P_protein_ttrR', 'ku'): 11.01321586, #this is the polymerase unbinding rate
('transcription_mm', 'ktx'): .01, #These are the parameters for transcription leak
('transcription_mm','kb'):100, #default binding rate!
#('ttrS-P',"kdephos"):1000, #auto dephosphorylation
('ligand_tt_protein_ttrS-P','kdephos'):1000, #dephosphorylation rate of ttrS
#('ttrS-P',"kphos"):10, #ttrs phosphorylates itself automatically. This is a proxy for sensing
('ligand_tt_protein_ttrS-P','kphos'):10, #equivalent to kphos, from before
#('ttrS-ttrR','ktransfer'):5, #rate of transfer from ttrS to ttrR
('ligand_tt_protein_ttrS-ttrR', 'ktransfer'):5, #this is the new transfer rate
('ttrR-P',"kdephos"):1000, #auto dephosphorylation
"kb":100, "ku":10, "ktx":.05, "ktl":.2, "kdeg":3, 'cooperativity':2, #some default parameters
('phrpL_hrpR', 'cooperativity'): 1.0, #default is 2, but these proteins bind at cooperativity = 1
('phrpL_hrpS', 'cooperativity'): 1.0,
('phrpL_hrpR_hrpS_RNAP', 'ktx'): 0.17025, #this is the transcription rate
('phrpL_hrpR_hrpS_RNAP', 'ku'): 11.01321586, #this is the polymerase unbinding rate
('translation_mm', 'B0030', 'ku'): 888, #Unbinding
('translation_mm', 'B0030', 'ktl'): .15, #Translation Rate
('transcription_mm', 'ktx'): .01, #These are the parameters for transcription leak
('transcription_mm', 'ku'): 100,
('phrpL_hrpR','ku'):100,
('phrpL_hrpS','ku'):100,
('phrpL_hrpS_RNAP', 'ku'): 1000.,#we'll say that leak has 100x the unbinding rate
('phrpL_hrpR_RNAP', 'ku'): 1000,
('phrpL','ku'):100000,
('binding', 'pLac_lacI','ktx'):0.001,
('plac_lacI','kb'):10,
('plac_lacI','ku'):0.001,
('binding', 'plac_hrpS_2x_lacI','kb'):8888,
("simple_transcription", 'plac','ktx'):0.1,
#Make a lot of repressor
('translation_mm', 'rbsL', 'ku'): 0.0001, #Unbinding
('translation_mm', 'rbsL', 'ktl'): 2, #Translation Rate
('plac_hrpS','kb'):100,
('plac_hrpS','ku'):0.0001
}
#P_rbsS_ttrS
constit_ttrS = DNAassembly(name = "constit_ttrS", promoter = "P", rbs = "rbsS", protein = "ttrS")
#P_rbsR_ttrR
constit_ttrR = DNAassembly(name = "constit_ttrR", promoter = "P", rbs = "rbsR", protein = "ttrR")
#P_arl_lacI
constit_lacI = DNAassembly(name = "constit_lacI", promoter = "P", rbs = "rbsL", protein = "lacI")
#P_arl_lacI
constit_mScarlet = DNAassembly(name = "constit_mScarlet", promoter = "P", rbs = "B0034", protein = "mScarlet")
# plac_hrpS
IPTG = Species("IPTG", material_type='inducer')
#IPTG = Protein("IPTG")
lac_p = LacIPromoter2(name="pLac", IPTG = "IPTG", lacI = "lacI", leak=False)
plac_hrpS = DNAassembly(name = "plac_hrpS", promoter=lac_p, rbs="rbsH",protein="hrpS")
# pttr_hrpR
tt = Species("tt",material_type="ligand")
ttrSbound = ChemicalComplex([Protein("ttrS"),tt])
ttrS3 = PhosphoTransferase(ttrSbound.species,targets=["ttrR"])
ttrR = PhosphoProtein("ttrR")
# Need to change from UTR1
pttr_hrpR = DNAassembly(name = "pttr_hrpR",promoter=RegulatedPromoter("pttr",[ttrR.get_phosphorylated_species()],leak=False),rbs="UTR1",protein="hrpR")
hrpL_gfp = DNAassembly("hrpL_gfp",promoter=CombinatorialPromoter("phrpL",["hrpR","hrpS"], tx_capable_list = [['hrpS', 'hrpR'], ['hrpS'], ['hrpR']], leak=False),rbs="B0030",protein="GFP")
extract_1_TXTL = TxTlExtract(name = "e coli", components = [pttr_hrpR,ttrSbound,ttrS3,ttrR,constit_ttrS,constit_ttrR,
constit_lacI, plac_hrpS, hrpL_gfp, constit_mScarlet],parameters=parameters)
CRN_extract_1 = extract_1_TXTL.compile_crn()
param_printer2(x0, parameters)
x0 = {"dna_pttr_hrpR":1,
"dna_constit_ttrS":1,
"dna_constit_ttrR":1,
"dna_constit_lacI":1,
"dna_plac_hrpS":1,
"dna_hrpL_gfp":1,
"dna_constit_mScarlet":5,
"phosphate_P":100,
"protein_RNAP":15,
"protein_RNAase":30,
"protein_Ribo":120,
"protein_lacI": 5}
IPTG_um_list = np.logspace(-3, 2, 6).round(decimals=6)
tt_um_list = np.logspace(-3, 2, 6).round(decimals=6)
fig, heat_map = simulate_heatmap('ligand_tt',tt_um_list,'inducer_IPTG',IPTG_um_list,x0=x0,normalize=True,title='Compiled Full AND gate \n lacI Spike', CRN=CRN_extract_1)
heat_map.set(xlabel='IPTG (uM)', ylabel='Tetrathionate (uM)')
heat_map.invert_yaxis()
plt.yticks(fontsize=18, fontstyle='normal', rotation = 0)
fig.savefig("./figures/fig6/fig6d.png")
plt.show()
x0 = {"dna_pttr_hrpR":1,
"dna_constit_ttrS":1,
"dna_constit_ttrR":1,
"dna_constit_lacI":1,
"dna_plac_hrpS":1,
"dna_hrpL_gfp":1,
"dna_constit_mScarlet":5,
"phosphate_P":100,
"protein_RNAP":15,
"protein_RNAase":30,
"protein_Ribo":120,
"protein_lacI": 0}
IPTG_um_list = np.logspace(-3, 2, 6).round(decimals=6)
tt_um_list = np.logspace(-3, 2, 6).round(decimals=6)
fig, heat_map = simulate_heatmap('ligand_tt',tt_um_list,'inducer_IPTG',IPTG_um_list,x0=x0,normalize=True,title='Compiled Full AND gate \n No lacI Spike', CRN=CRN_extract_1)
heat_map.set(xlabel='IPTG (uM)', ylabel='Tetrathionate (uM)')
heat_map.invert_yaxis()
plt.yticks(fontsize=18, fontstyle='normal', rotation = 0)
fig.savefig("./figures/fig6/fig6b.png")
plt.show()
timepoints = np.linspace(0, 200, 1000)
x0 = {"dna_pttr_hrpR":5,
"dna_constit_ttrS":1,
"dna_constit_ttrR":1,
# Cole origin
"dna_constit_lacI":40,
"dna_plac_hrpS":1,
"dna_hrpL_gfp":1,
"dna_constit_mScarlet":5,
"phosphate_P":100,
"protein_RNAP":10,
"protein_Ribo":90.,
"ligand_tt": 1,
"inducer_IPTG": 0.,
"protein_lacI": 0,}
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
results_melted = pd.melt(Re1, var_name='species', id_vars=['time'])
results_melted.head()
Re1 = CRN_extract_1.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)
p = bokeh.plotting.figure(title = 'Early Circuit Dynamics, pLac', x_axis_label = 'Time (min)', y_axis_label = 'Simulated Concentration (uM)', height = 2000, width = 4300)
p.line(timepoints, Re1['complex_dna_plac_hrpS_protein_RNAP'], legend_label = 'Active pLac', color = 'Green', line_width = 14)
p.line(timepoints, Re1['rna_constit_lacI'], legend_label = 'LacI RNA', color = 'orange', line_width = 14)
p.line(timepoints, Re1['protein_hrpS'], legend_label = 'hrpS protein', color='lightsteelblue', line_width = 14)
p.line(timepoints, Re1['complex_dna_hrpL_gfp_protein_hrpR_protein_hrpS_protein_RNAP'], legend_label = 'hrpS, hrpR both bound to pHRPL', color='darkorchid', line_width = 14)
p.line(timepoints, Re1['rna_hrpL_gfp'], legend_label = 'GFP RNA', color='navy', line_width = 14)
p.add_layout(p.legend[0], 'right')
p.legend.title = 'Species'
p.legend.glyph_width = 100
p.legend.padding = 100
p.legend.title_text_font_size = '80pt'
p.legend.label_standoff = 50
p.legend.glyph_height = 100
p.xaxis.major_tick_line_width = 10
p.yaxis.major_tick_line_width = 10
p.xaxis.minor_tick_line_width = 10
p.yaxis.minor_tick_line_width = 10
p.axis.minor_tick_out = 20
p.axis.minor_tick_out = 20
p.axis.major_tick_out = 40
p.axis.major_tick_out = 40
p.yaxis.axis_label_standoff = 50
p.xaxis.axis_label_standoff = 50
p.xaxis.axis_line_width = 8
p.yaxis.axis_line_width = 8
p.yaxis.major_label_standoff = 30
p.xaxis.major_label_standoff = 30
p.ygrid.grid_line_width = 10
p.xgrid.grid_line_width = 10
p.axis.axis_label_text_font_size = "80pt"
p.legend.label_text_font_size = "80pt"
p.title.text_font_size = '100pt'
p.yaxis.major_label_text_font_size = "80pt"
p.xaxis.major_label_text_font_size = "80pt"
#p.y_range = Range1d(0, 0.9)
bokeh.io.show(p)
%load_ext watermark
%watermark -v -p numpy,pandas,scipy,bioscrape,biocrnpyler,bokeh,jupyterlab,holoviews